home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / LSODE.cc < prev    next >
C/C++ Source or Header  |  1996-07-24  |  9KB  |  456 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <cfloat>
  32. #include <cmath>
  33.  
  34. #include <iostream.h>
  35.  
  36. #include "LSODE.h"
  37. #include "f77-fcn.h"
  38. #include "lo-error.h"
  39.  
  40. extern "C"
  41. {
  42.   int F77_FCN (lsode, LSODE) (int (*)(const int&, const double&,
  43.                       double*, double*, int&),
  44.                   int&, double*, double&, double&, int&,
  45.                   double&, double&, int&, int&, int&,
  46.                   double*, int&, int*, int&,
  47.                   int (*)(const int&, const double&,
  48.                       double*, const int&, const int&,
  49.                       double*, const int&),
  50.                   int&);
  51. }
  52.  
  53. static ODEFunc::ODERHSFunc user_fun;
  54. static ODEFunc::ODEJacFunc user_jac;
  55. static ColumnVector *tmp_x;
  56.  
  57. LSODE::LSODE (void) : ODE (), LSODE_options ()
  58. {
  59.   n = size ();
  60.  
  61.   stop_time_set = 0;
  62.   stop_time = 0.0;
  63.  
  64.   integration_error = 0;
  65.   restart = 1;
  66.  
  67.   istate = 1;
  68.   itol = 1;
  69.   itask = 1;
  70.   iopt = 0;
  71.  
  72.   liw = 20 + n;
  73.   lrw = 22 + n * (9 + n);
  74.  
  75.   sanity_checked = 0;
  76. }
  77.  
  78. LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f)
  79.   : ODE (state, time, f), LSODE_options ()
  80. {
  81.   n = size ();
  82.  
  83.   stop_time_set = 0;
  84.   stop_time = 0.0;
  85.  
  86.   integration_error = 0;
  87.   restart = 1;
  88.  
  89.   istate = 1;
  90.   itol = 1;
  91.   itask = 1;
  92.   iopt = 0;
  93.  
  94.   liw = 20 + n;
  95.   lrw = 22 + n * (9 + n);
  96.  
  97.   sanity_checked = 0;
  98. }
  99.  
  100. void
  101. LSODE::force_restart (void)
  102. {
  103.   restart = 1;
  104. }
  105.  
  106. void
  107. LSODE::set_stop_time (double time)
  108. {
  109.   stop_time_set = 1;
  110.   stop_time = time;
  111. }
  112.  
  113. void
  114. LSODE::clear_stop_time (void)
  115. {
  116.   stop_time_set = 0;
  117. }
  118.  
  119. int
  120. lsode_f (const int& neq, const double& time, double *,
  121.      double *deriv, int& ierr) 
  122. {
  123.   ColumnVector tmp_deriv;
  124.  
  125.   // NOTE: this won't work if LSODE passes copies of the state vector.
  126.   //       In that case we have to create a temporary vector object
  127.   //       and copy.
  128.  
  129.   tmp_deriv = (*user_fun) (*tmp_x, time);
  130.  
  131.   if (tmp_deriv.length () == 0)
  132.     ierr = -1;
  133.   else
  134.     {
  135.       for (int i = 0; i < neq; i++)
  136.     deriv [i] = tmp_deriv.elem (i);
  137.     }
  138.  
  139.   return 0;
  140. }
  141.  
  142. int
  143. lsode_j (const int& neq, const double& time, double *,
  144.      const int&, const int&, double *pd, const int& nrowpd)
  145. {
  146.   Matrix tmp_jac (neq, neq);
  147.  
  148.   // NOTE: this won't work if LSODE passes copies of the state vector.
  149.   //       In that case we have to create a temporary vector object
  150.   //       and copy.
  151.  
  152.   tmp_jac = (*user_jac) (*tmp_x, time);
  153.  
  154.   for (int j = 0; j < neq; j++)
  155.     for (int i = 0; i < neq; i++)
  156.       pd [nrowpd * j + i] = tmp_jac (i, j);
  157.  
  158.   return 0;
  159. }
  160.  
  161. ColumnVector
  162. LSODE::do_integrate (double tout)
  163. {
  164.   ColumnVector retval;
  165.  
  166.   if (restart)
  167.     {
  168.       restart = 0;
  169.       istate = 1;
  170.     }
  171.  
  172.   if (iwork.length () != liw)
  173.     {
  174.       iwork.resize (liw);
  175.  
  176.       for (int i = 4; i < 9; i++)
  177.     iwork.elem (i) = 0;
  178.     }
  179.  
  180.   if (rwork.length () != lrw)
  181.     {
  182.       rwork.resize (lrw);
  183.  
  184.       for (int i = 4; i < 9; i++)
  185.     rwork.elem (i) = 0;
  186.     }
  187.  
  188.   if (jac)
  189.     method_flag = 21;
  190.   else
  191.     method_flag = 22;
  192.  
  193.   integration_error = 0;
  194.  
  195.   double *xp = x.fortran_vec ();
  196.  
  197.   // NOTE: this won't work if LSODE passes copies of the state vector.
  198.   //       In that case we have to create a temporary vector object
  199.   //       and copy.
  200.  
  201.   tmp_x = &x;
  202.   user_fun = function ();
  203.   user_jac = jacobian_function ();
  204.  
  205.   if (! sanity_checked)
  206.     {
  207.       ColumnVector xdot = (*user_fun) (x, t);
  208.  
  209.       if (x.length () != xdot.length ())
  210.     {
  211.       (*current_liboctave_error_handler)
  212.         ("lsode: inconsistent sizes for state and derivative vectors");
  213.  
  214.       integration_error = 1;
  215.       return retval;
  216.     }
  217.  
  218.       sanity_checked = 1;
  219.     }
  220.  
  221.   // Try 5000 steps before giving up.
  222.  
  223.   iwork.elem (5) = 5000;
  224.  
  225.   if (stop_time_set)
  226.     {
  227.       itask = 4;
  228.       rwork.elem (0) = stop_time;
  229.     }
  230.   else
  231.     {
  232.       itask = 1;
  233.     }
  234.  
  235.   double abs_tol = absolute_tolerance ();
  236.   double rel_tol = relative_tolerance ();
  237.  
  238.   rwork.elem (4) = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0;
  239.   rwork.elem (5) = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0;
  240.   rwork.elem (6) = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0;
  241.  
  242.   int *piwork = iwork.fortran_vec ();
  243.   double *prwork = rwork.fortran_vec ();
  244.  
  245.   working_too_hard = 0;
  246.  
  247.  again:
  248.  
  249.   F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol,
  250.                abs_tol, itask, istate, iopt, prwork, lrw,
  251.                piwork, liw, lsode_j, method_flag));
  252.  
  253.   if (f77_exception_encountered)
  254.     (*current_liboctave_error_handler) ("unrecoverable error in lsode");
  255.   else
  256.     {
  257.       switch (istate)
  258.     {
  259.     case -13: // Return requested in user-supplied function.
  260.     case -6:  // error weight became zero during problem. (solution
  261.               // component i vanished, and atol or atol(i) = 0.)
  262.     case -5:  // repeated convergence failures (perhaps bad jacobian
  263.               // supplied or wrong choice of mf or tolerances).
  264.     case -4:  // repeated error test failures (check all inputs).
  265.     case -3:  // illegal input detected (see printed message).
  266.     case -2:  // excess accuracy requested (tolerances too small).
  267.       integration_error = 1;
  268.       break;
  269.  
  270.     case -1:  // excess work done on this call (perhaps wrong mf).
  271.       working_too_hard++;
  272.       if (working_too_hard > 20)
  273.         {
  274.           (*current_liboctave_error_handler)
  275.         ("giving up after more than %d steps attempted in lsode",
  276.          iwork.elem (5) * 20);
  277.           integration_error = 1;
  278.         }
  279.       else
  280.         {
  281.           istate = 2;
  282.           goto again;
  283.         }
  284.       break;
  285.  
  286.     case 2:  // lsode was successful
  287.       retval = x;
  288.       t = tout;
  289.       break;
  290.       
  291.     default: // Error?
  292.       (*current_liboctave_error_handler)
  293.         ("unrecognized value of istate returned from lsode");
  294.       break;
  295.     }
  296.     }
  297.  
  298.   return retval;
  299. }
  300.  
  301. #if 0
  302. void
  303. LSODE::integrate (int nsteps, double tstep, ostream& s)
  304. {
  305.   int time_to_quit = 0;
  306.   double tout = t;
  307.  
  308.   s << t << " " << x << "\n";
  309.  
  310.   for (int i = 0; i < nsteps; i++)
  311.     {
  312.       tout += tstep;
  313.       if (stop_time_set && tout > stop_time)
  314.     {
  315.       tout = stop_time;
  316.       time_to_quit = 1;
  317.     }
  318.  
  319.       x = integrate (tout);
  320.  
  321.       s << t << " " << x << "\n";
  322.  
  323.       if (time_to_quit)
  324.     return;
  325.     }
  326. }
  327. #endif
  328.  
  329. Matrix
  330. LSODE::do_integrate (const ColumnVector& tout)
  331. {
  332.   Matrix retval;
  333.   int n_out = tout.capacity ();
  334.  
  335.   if (n_out > 0 && n > 0)
  336.     {
  337.       retval.resize (n_out, n);
  338.  
  339.       for (int i = 0; i < n; i++)
  340.     retval.elem (0, i) = x.elem (i);
  341.  
  342.       for (int j = 1; j < n_out; j++)
  343.     {
  344.       ColumnVector x_next = do_integrate (tout.elem (j));
  345.  
  346.       if (integration_error)
  347.         return retval;
  348.  
  349.       for (int i = 0; i < n; i++)
  350.         retval.elem (j, i) = x_next.elem (i);
  351.     }
  352.     }
  353.  
  354.   return retval;
  355. }
  356.  
  357. Matrix
  358. LSODE::integrate (const ColumnVector& tout, const ColumnVector& tcrit)
  359. {
  360.   Matrix retval;
  361.   int n_out = tout.capacity ();
  362.  
  363.   if (n_out > 0 && n > 0)
  364.     {
  365.       retval.resize (n_out, n);
  366.  
  367.       for (int i = 0; i < n; i++)
  368.     retval.elem (0, i) = x.elem (i);
  369.  
  370.       int n_crit = tcrit.capacity ();
  371.  
  372.       if (n_crit > 0)
  373.     {
  374.       int i_crit = 0;
  375.       int i_out = 1;
  376.       double next_crit = tcrit.elem (0);
  377.       double next_out;
  378.       while (i_out < n_out)
  379.         {
  380.           int do_restart = 0;
  381.  
  382.           next_out = tout.elem (i_out);
  383.           if (i_crit < n_crit)
  384.         next_crit = tcrit.elem (i_crit);
  385.  
  386.           int save_output;
  387.           double t_out;
  388.  
  389.           if (next_crit == next_out)
  390.         {
  391.           set_stop_time (next_crit);
  392.           t_out = next_out;
  393.           save_output = 1;
  394.           i_out++;
  395.           i_crit++;
  396.           do_restart = 1;
  397.         }
  398.           else if (next_crit < next_out)
  399.         {
  400.           if (i_crit < n_crit)
  401.             {
  402.               set_stop_time (next_crit);
  403.               t_out = next_crit;
  404.               save_output = 0;
  405.               i_crit++;
  406.               do_restart = 1;
  407.             }
  408.           else
  409.             {
  410.               clear_stop_time ();
  411.               t_out = next_out;
  412.               save_output = 1;
  413.               i_out++;
  414.             }
  415.         }
  416.           else
  417.         {
  418.           set_stop_time (next_crit);
  419.           t_out = next_out;
  420.           save_output = 1;
  421.           i_out++;
  422.         }
  423.  
  424.           ColumnVector x_next = do_integrate (t_out);
  425.  
  426.           if (integration_error)
  427.         return retval;
  428.  
  429.           if (save_output)
  430.         {
  431.           for (int i = 0; i < n; i++)
  432.             retval.elem (i_out-1, i) = x_next.elem (i);
  433.         }
  434.  
  435.           if (do_restart)
  436.         force_restart ();
  437.         }
  438.     }
  439.       else
  440.     {
  441.       retval = do_integrate (tout);
  442.  
  443.       if (integration_error)
  444.         return retval;
  445.     }
  446.     }
  447.  
  448.   return retval;
  449. }
  450.  
  451. /*
  452. ;;; Local Variables: ***
  453. ;;; mode: C++ ***
  454. ;;; End: ***
  455. */
  456.